# TODO: Add comment
# 
# Author: broken
###############################################################################
#import library
library(base)
library(raster)
library(gstat)
library(automap)
library(stringr)


#path
#regressPm_file="/home/phamha/Desktop/npp.tif"
regressPm_file = "C:/npp_2.tif"

	
	#PM values
	pmRaster=raster(regressPm_file)
	pm=values(pmRaster)
	corxy=coordinates(pmRaster)
	x=corxy[,'x']
	y=corxy[,'y']
	

	totalCell=length(pmRaster)
	cell = c(1:totalCell)

	table=data.frame(cell,x,y,pm)
	newTable=table
	trainTable=subset(table,!is.na(pm)&pm!=-9999)
	testTable=subset(table,is.na(pm))

	auto_trainTable = trainTable
	coordinates(auto_trainTable) =~ x+y
	
	auto_variogram = autofitVariogram(pm~1,auto_trainTable)
	auto_model = auto_variogram$var_model$model[2]
	auto_sill = auto_variogram$var_model$psill[2]
	auto_range = auto_variogram$var_model$range[2]
	
	#caculate variogram
	empiVario = variogram(pm~1,locations=~x+y,data=trainTable)
	
	sphModel=vgm(psill=auto_sill,model=auto_model,nugget=0,range=auto_range)		
	sphFit=fit.variogram(empiVario,sphModel)
	#plot(empiVario,model=sphFit)
	#sph fit
	#sill=min(empiVario$gamma)
	#sphModel=vgm(psill=sill,model="Sph",nugget=0,range=min(empiVario$dist))
	#sphModel=vgm(model="Sph",nugget=0,range=1)		
	#sphFit=fit.variogram(empiVario,sphModel)
	
	universal_result=krige(id="pm",formula=pm~x+y,data=trainTable,newdata=testTable,model=sphFit,locations=~x+y)
	

	#edit tiff
	newTable$pm[is.na(newTable$pm)] = universal_result[,3]
	universalPMRaster=pmRaster
	universalPMRaster[1:totalCell]=newTable$pm
	
	## #edit error tiff
	## errorPMRaster=pmRaster
	## errorPMValue=universal_result[,4]
	## errorPMRaster[1:totalCell]=errorPMValue
	
	#save uk result to tiff
	uk_file = str_replace(regressPm_file,".tif","_uk.tif")
	writeRaster(universalPMRaster,filename=uk_file,format="GTiff")
	
	#save uk error to tiff
	#error_file = str_replace(regressPm_file,".tif","_error.tif")
	#writeRaster(errorPMRaster,filename=error_file,format="GTiff")
	print(uk_file)
	







